Strongly Coupled Plasmon Polaritons in Gold and Epsilon-Near-Zero Bifilms

Epsilon-near-zero (ENZ) polaritons in a thin transparent conducting-oxide film exhibit a significant electric field enhancement and localization within the film at frequencies close to their plasma frequency, but do not propagate. Meanwhile, plasmon polariton modes in thin metallic films can propagate for several microns, but are more loosely confined in the metal. Here, we propose a strongly coupled bilayered structure of a thin gold film on a thin indium tin oxide (ITO) film that supports hybrid polariton modes. We experimentally characterize the dispersion of these modes and show that they have propagation lengths of 4–8 μm while retaining mode confinement greater than that of the polariton in gold films by nearly an order of magnitude. We study the tunability of this coupling strength by varying the thickness of the ITO film and show that ultrastrong coupling is possible at certain thicknesses. The unusual linear and nonlinear optical properties of ITO at ENZ frequencies make these bifilms useful for the active tuning of strong coupling, ultrafast switching, and enhanced nonlinear interactions at near-infrared frequencies.

sample assembly is mounted on a rotation mount to vary the incidence angle on the sample.
Another lens L2 images the sample onto the entrance pupil of a microscope objective O2, which couples the reflected light to a multi-mode fiber coupled spectrometer. We use an optical spectrum analyzer OSA (Agilent 86142A) to record the spectra from 600 nm to 1600 nm. We use an InGaAs spectrometer (customized SM304 from Spectral Products) to record the spectra from 1600 nm to 2300 nm.
We take into account the shape of the blackbody spectrum of the lamp, the spectral response of the prisms and the oil, and the optical elements used in the setup by normalizing the TM-polarized reflectance spectrum from the sample at a particular incidence angle to the TE-polarized reflectance spectrum at the same angle. Before the spectral normalization, the detector noise in each spectrum is smoothed over first by using a Savtizky-Golay filter, and then by applying the wavelet transform.

S3. Reflectance maps of the standalone constituent modes
The reflectance (R T M /R T E ) maps of the standalone gold sample are shown in the top panels in Figure S3, while the reflectance maps of the standalone 23 nm thick ITO (permittivity spectrum shown in fig S1( Figure S4: Reflectance map R TM /R TE obtained from TMM simulations in the un-normalized wavevector and frequency space of (a) a bifilm with 35-nm-thick gold and 23-nm-thick ITO with the same permittivity as in bifilm A, and of (b) bifilm A, and (c) bifilm C with reduced losses in the ITO film. The spectra of ITO used in (b) and (c) are shown in Figures S1(a) and S1(c), respectively but with Im[ ITO ] reduced by a factor of 10. Asymmetric lorentzian fits (red, dashed) to the extinction spectrum ((1 − R TM /R TE ), blue, solid) of bifilm A for the (c) lower, and the (d) upper polariton at a certain Re[k x ] (shown in the inset of each plot). The 95% confidence intervals of the fit parameters [ω 0 /(2π, γ)] are stated in the inset.
The critical coupling condition for our bifilm structure at each frequency ν, or when the coupling losses are balanced by the absorption losses, 1 are given by the solutions to the characteristic equation in the complex wavevector k x and real frequency ν space for which the Im[k x ] is minimized. These solutions depend on the geometrical parameters of the bifilm. The same set of parameters do not satisfy the critical coupling criterion for both the upper and the lower polariton simultaneously due to large polariton band gap. Figure S4(a) shows the reflectance map of a bifilm with a 35-nm-thick gold film and 23-nm-thick ITO film obtained from TMM simulations in the frequency ν and un-normalized wavevector k x space. The permittivity of the ITO layer is taken to be the same as bifilm A. Comparing S4 Figure S4(a) with the reflectance map of bifilm A shown in Figure 3(a) of the main text, we note that the choice of a thinner gold film in the former case leads to a better coupling efficiency to the lower polariton. However, the thinner gold film also leads to larger radiative damping of the upper polariton, which is evident in its broader spectral linewidth. Hence, the choice of a 50-nm-thick gold layer for our bifilms leads to efficient (inefficient) coupling to the upper (lower) polariton.
Large absorption losses within the ITO film also contribute to the smaller amplitude of the resonance dip for the lower polariton compared to the upper polariton. Figure S4 Figure S2). The reduced losses in the ITO film lead to much sharper dips for both polariton branches when compared with the reflectance map calculated with the full permittivity of the ITO film ITO (Figure 3(a) in the main text).
Additionally, the reduced losses in ITO result in a well-defined lower polariton branch that continues to exist until it approaches the prism light line. Figure S4(c) shows the reflectance map of bifilm C with similarly reduced losses in the ITO layer, and we see a well defined lower polariton branch that is pushed to smaller frequencies than for the thinner ITO in where Here A is a fitting parameter that determines the peak of the extinction spectrum, and a is the asymmetry parameter. When a is zero, we recover the normal symmetric Lorentzian function with a full width at half maximum linewidth of 2γ 0 . The fit parameters B and C account for the frequency dependent linear distortion of the extinction spectra due to the dispersion of the prism.
From the dispersion lines of the SPP modeω SPP (k x ), and the upper and the lower hybrid polaritonsω U,L (k x ) calculated from their respective reflectance maps, we estimate the coupling strength g R of the bifilm system by fitting the eigenvalues of the Hopfield-Bogliubov interaction Hamiltonian matrix, 3 shown in equation (3), to the upper and the lower polariton We assume that the ENZ mode has a flat dispersion line given bỹ where γ ITO is the damping in the Drude permittivity model of ITO, and the resonance S6 frequency ω 0,ENZ is close to the ENZ frequency of ITO given by where ω P is the plasma frequency, and ∞ is the asymptotic value of permittivity for frequencies much larger than the ENZ frequency in the Drude permittivity model of ITO. We take ω 0,ENZ to be an adjustable parameter in our fit, and not equal to ω ENZ,ITO as this assumption only holds true for an infinitesimally small thickness of the ITO film. 4 As the ITO thickness increases, the symmetric polariton mode within the ITO film transitions from an ENZ mode, which is asymptotically pinned to ω ENZ,ITO , to a long-range surface plasmon polariton (LR-SPP) mode, whose dispersion line asymptotically approaches lower frequencies than ω ENZ,ITO for large wavevectors.
We begin with ω 0,ENZ = ω ENZ,ITO , and use the nonlinear least squares method to individually fit the dispersion lines of the upper and the lower polaritons to their analytical dispersion relations given in equation (3). If the difference between the two values of g R that we obtain from the upper and the lower polariton fits is larger than 10%, we repeat the curve fit albeit with a slightly reduced value of ω 0,ENZ until the difference in g R obtained from both the upper and the lower polariton fits is minimized. The average of the two g R values obtained after this optimization procedure is the final calculated value of the coupling strength of the hybrid polaritons, and the values themselves form the confidence interval of the fit.
We note that as the thickness of the ITO film d ITO is increased, the assumption of a flat dispersion line for the ENZ mode at all wavevectors no longer holds true as the mode becomes more LR-SPP like in nature, even with the optimization of its resonance frequency.
This effect manifests in relatively larger confidence intervals of g R in the fit (shown by the shaded blue region in Figure 3(d) of the main text) for d ITO > 65 nm than for bifilms with thinner ITO. Additionally, g R itself also increases as d ITO is increased due to an enhancement S7 of the oscillator strength of the ENZ mode. Consequently, as discussed in the main text, g R becomes larger than 0.1ω 0,ENZ , where ω 0,ENZ is the avoided crossing frequency, and the SPP and the ENZ modes in the bifilm are ultra-strongly coupled for d ITO larger than 30 nm.
In this ultra-strong coupling regime, the simple analytical (Hopfield) model for the hybrid polaritons based on two coupled harmonic oscillators is not accurate due to several underlying assumptions made by the model. For instance, the self-interaction term of the SPP and the ENZ modes, which is proportional to the square of their respective polarizations, is neglected in the interaction Hamiltonian. 5 In the complete interaction Hamiltonian written in the dipole gauge, this self-interaction term is responsible for the renormalization of the respective uncoupled mode frequencies.
In the quantum-mechanical picture, this term corresponds to the A 2 term of the interaction Hamiltonian, which also includes the counter-rotating terms that are neglected when making the rotating-wave approximation (RWA). This selfinteraction term leads to deviations from the curves obtained from the Hopfield model. For intersubband polaritons formed by ultra-strong coupling between a microcavity mode and the bound states in multiple quantum wells (which form a collective Berreman mode for S8 a large enough number of quantum wells), it has been previously demonstrated that the lowenergy polariton branch asymptotically approaches a smaller energy than the high-energy polariton close to the avoided-crossing. 6 As a consequence, the high reflectivity between the two energy asymptotes opens up a so-called "polaritonic gap" in the dispersion line of the intersubband polaritons. We observe a similar effect in the dispersion lines of our bifilms that support ultra-strongly coupled SPP and ENZ modes, that is for bifilms with d ITO > 30 nm. As shown in Figure S5, the frequency asymptote of the lower polariton (blue, solid) is smaller than the frequency at that specific transverse wavevector in the fitted dispersion where l = {s, i, a, p} is the index for the substrate, ITO, gold and prism layers, respectively, in this multi-layered structure; k x is the transverse wavevector; a l and b l are the coefficients of the forward and backward propagating solutions within the layer that are determined by the field continuity relations, and k zl is the longitudinal wavevector given by with l being the permittivity of the layer l. The continuity of E xl and D zl (= l E zl ) at each interface leads to a set of 8 linear homogeneous equations for the 8 field coefficients In matrix form, the set of equations can be S10 written as where L is given by   Figure S1(a)) on glass, respectively at their point of degeneracy, or where their respective dispersion lines cross. As mentioned previously, the SPP mode is mostly confined at the gold-substrate interface, with a very small longitudinal field amplitude E z in the gold layer. The ENZ mode, on the other hand, is tightly confined within the ITO film, where E z is enhanced almost by a factor of three than at the ITO-substrate interface, and is largely constant within the film.    Figure S7(b), respectively at their point of degeneracy. We see that for this thicker ITO film, the "ENZ" mode starts to resemble the SPP mode so that the fields S14 are now confined along the ITO-substrate interface, and E z is no longer strongly enhanced within the ITO film. For the bifilm comprising the 50 nm thick gold on 65 nm thick ITO, the electric field of the hybrid polaritons within the ITO film starts to decouple between the two ITO interfaces, as seen in the electric field profiles of the upper and the lower polariton close to the avoided crossing region shown in Figures S8(c) and S8(d), respectively. Both For even thicker ITO films, the "ENZ" mode becomes even more SPP-like in nature. As shown in Figure S9 To calculate the propagation lengths of the hybrid polaritons, and the SPP mode from S16 their reflectance map -both simulated and measured -we perform a nonlinear least-squares curve fit of asymmetric Lorentzian functions of the following form to the extinction (1 − R TM /R T E ) wavevector scan at each frequency where The fitting function is similar to the one assumed for the extracting the decay rates γ 0 from the reflectance maps. The fit parameter k x0 calculated from the wavevector scans at various ω produces the dispersion line, albeit with back-bending present where relevant. 8,9 On the other hand, the propagation length L prop of the polariton is given by To calculate L prop of the ENZ mode, we first extract the linewidth γ of the spectral dip from the reflectance maps shown in Figures S2. We then use the following relation from the S17 supplement in Ref. 10 to calculate the propagation length  respectively). However, the slope of the confinement curve is steeper for bifilm B than for bifilm A, and becomes as large as 0.4 for bifilm B away from the polariton band edge, which is more than thrice the confinement at that particular frequency (350 THz) for the upper polariton in bifilm A. The mode confinement of the upper polariton in bifilm C is even larger than the upper polariton in bifilm B throughout spectral range with a maximum value close to 0.5. This trend is in agreement with the transition of the hybrid polaritons to interface polaritons for larger thicknesses of the ITO layer discussed previously in section S6. At the "bluer" frequencies, where ITO behaves as a dielectric, the (upper polariton) becomes S19 more confined along the gold-ITO interface for both bifilms B and C. This confinement is even larger for bifilm C where the ENZ mode in the 100 nm-thick-ITO layer is even more LR-SPP-like than the ENZ mode in the 65 nm-thick-ITO layer in bifilm B.
On the other hand, the longitudinal field enhancement of the upper polariton in bifilm A is larger than both bifilms B and C throughout the spectral range under consideration. The maximum field enhancement for the upper polariton in bifilm A ≈ 32× at the band edge of bifilm A, whereas it is ≈ 20× and ≈ 17× at the band edges of bifilms B and C, respectively.
A larger field enhancement would be more desirable for enhancing the nonlinear optical response for certain applications. The larger mode confinements of the upper polaritons in bifilms B and C in comparison to the polariton in bifilm A also reflects in their smaller propagation lengths and larger damping. The simulated as well as the measured propagation lengths for the upper polaritons in both bifilms B and C is between 2-4 µm, whereas the corresponding lengths for bifilm A is between 4-8 µm. On the other hand, the simulated and the measured damping rates of the upper polaritons in both bifilms B and C is larger than 0.1γ ITO , while the damping rate for the upper polariton in bifilm A is lower than 0.1γ ITO throughout.